%% fig 5b; place cell recruitment   

clear
load('/media/Big1/15tracks/recruitment/recruitment_newrand.mat','out');
savepath='/media/Big1/15tracks/recruitment/';

load('/media/Big1/15tracks/recruitment/recruitment_eachpeakRate_newrand.mat');
%
 yall=cellfun(@(x) x{1}(:),{out.parti_y},'UniformOutput',false);
 y=mean(cell2mat(yall),2);
%  y=y-y(1);
    x=out(1).trackx(:);
    %
    ntrack1celld1=cell2mat(cellfun(@(x) x(:),{out.ntr_r1d1},'UniformOutput',false)');
    ntrack1celld2=cell2mat(cellfun(@(x) x(:),{out.ntr_r1d2},'UniformOutput',false)');
    
    data1=ntrack1celld1;
    [gpParam1,poissParam1,dataHist1]=fit_gammapoiss_distribution(data1,1,{'#Track/cell','Cell'});
    data2=ntrack1celld2;
    [gpParam2,poissParam2,dataHist2]=fit_gammapoiss_distribution(data2,1,{'#Track/cell','Cell'});
   setaxisformalAll
   
% savepdf(gcf,'/media/Big1/15tracks/figures//PlFieldsGammapoissDistribut.pdf');
 
    figure;
    subplot(121)
    bar(dataHist1.edgex,dataHist1.count,1,'facealpha',0.5)
    hold on;
    bar(dataHist2.edgex,dataHist2.count,1,'facealpha',0.5)
   %
   peakRate=unique([out.peakRateThr]);
   for i=2%:length(peakRate) 
   data=cell2mat(cellfun(@(x) x(:),{out([out.peakRateThr]==peakRate(i)).ntr_r1d1},'UniformOutput',false)');
   
   distributionName={'NegativeBinomial','Poisson'};
    [disParam,dataHist]=fit_sample_distribution(data,distributionName,1,{'#Track/cell','Cell'});
    savepdf(gcf,fullfile(savepath,['DistributionOfNumberOfTrackPerCell',num2str(peakRate(i)),'Hz.pdf']));
   end
    setaxisformal
% savepdf(gcf,'/media/Big1/15tracks/figures/NumberOfTrackPerNeuronGammapoissDistribut.pdf');
%
    out1=out([out.peakRateThr]==2);
     yall=cellfun(@(x) x{1}(:),{out1.parti_y},'UniformOutput',false);
 y=mean(cell2mat(yall),2);
%  y=y-y(1);
    x=out1(1).trackx(:);
    xi=linspace(x(1),x(end),100);
    
    fun=@(pa,x) pa(3)*(x.^pa(1)./(x.^pa(1)+pa(2).^pa(1))); % Naka-Rushton function
    pa0=[.9,.9,101];
    palb=[-5,0.1,20];
    paub=[5,50,150];
    fitpa=lsqcurvefit(fun, pa0,x(2:end),y(2:end),palb,paub,optimoptions('lsqcurvefit',...
            'FunctionTolerance',1e-20,'StepTolerance',1e-30,'MaxIterations',1000,...
            'MaxFunctionEvaluations',1000,'OptimalityTolerance',1e-20)); 
       
%     fitpa=lsqcurvefit(fun,pa0,x,y,palb,paub,optimoptions('lsqcurvefit',...
%             'FunctionTolerance',1e-20,'StepTolerance',1e-20,'MaxIterations',1000,...
%             'MaxFunctionEvaluations',1000,'OptimalityTolerance',1e-20));
        
    yi=fun(fitpa,xi);
    
    targetPercent=95;
    xi100=x(1):1:1000;
    yi100=fun(fitpa,xi100);
    ind100=find(abs(yi100-targetPercent)==min(abs(yi100-targetPercent)));
    x100=xi100(ind100);
    xi100=xi100(max([1,ind100-20]):min([length(xi100),ind100+20]));
    yi100=yi100(max([1,ind100-20]):min([length(yi100),ind100+20]));
    
  
    figure('color','w','position',[472    84   626   615]);
    subplot(221)
    for i=1:length(out1)
       errorbar(out1(i).trackx,out1(i).parti_y{1},out1(i).parti_y{2},...
           'linewidth',1,'linestyle','-');hold on;
    end
    legend({out1.ratname},'location','best');
    xlabel('#Tracks');ylabel('%Cells');
    title('All rats')
    subplot(223)
    plot(x,y,'bo',xi,yi,'r-');
    xlabel('#Tracks');ylabel('%Cells');
    subplot(224)
    plot(xi100,yi100,'r-',xi100,targetPercent*ones(size(xi100)),'k-');
    title(x100); 
    subplot(222)
    for i=1:length(out1)
    plot(out1(i).dcell,'.-');hold on;
    end
    xlabel('Add 1 track');ylabel('Add % cells');
   
    setaxisformalAll
% % % %     savepdf(gcf,fullfile(savepath,'recruitCurve2Hz.pdf'));
%     savepdf(gcf,fullfile('/media/Big1/15tracks/figures/','PlacecellRecruitment_extrapolate.pdf'));
  % 
    eval_capacity_formula(y/100,x,{'linear','log','power','exp','hyperbolic'},1,{'#Track','%Cell'});
% % % % % % % % % % %      savepdf(gcf,fullfile(savepath,'compareModel.pdf'));
     setaxisformalAll
%      savepdf(gcf,fullfile('/media/Big1/15tracks/figures/','PlacecellRecruitment_compareModel.pdf'));
%% fig s5d compare with random sampling cells
clear
load('/media/Big1/15tracks/recruitment/recruitment_eachpeakRate_newrand.mat');
savepath='/media/Big1/15tracks/figures/';
out1=out([out.peakRateThr]==2);
 
ncell1tr=cell2mat(cellfun(@(x) sum(x),{out1.parti_matrix_2dir},'UniformOutput',false)');
[ncell1trm,ncell1trse]=mean_se(ncell1tr);
[ncell1trm0,ncell1trse0]=mean_se(ncell1tr(:));
p0=anova1(ncell1tr(:),reshape(repmat(1:15,5,1),5*15,1),'off');


[ncell1trm,ncell1trse]=mean_se(ncell1tr)

f3=figure('color','w');
shadedErrorBar(1:15,ncell1trm,ncell1trse,'b.-',0.6);
set(gca,'ylim',[0,40])
xlabel('Track');ylabel('# of neuron');
title(sprintf('p=%.2g, N=%.4g +- %.2g',p0,ncell1trm0,ncell1trse0))
setaxisformal

yall=cellfun(@(x) x{1}(:),{out1.parti_y},'UniformOutput',false);
y=mean(cell2mat(yall),2);
x=out(1).trackx(:);
 

yallncell=cellfun(@(x,y) x{1}(:)*y/100,{out1.parti_y},{out1.ncell_pyr},'UniformOutput',false);
yncellmat=cell2mat(yallncell);
yncell=mean(cell2mat(yallncell),2);

% ----by formula
randfun=@(x,m,N) (1-(1-m/N).^x)*N; 
yrandbyfun=cell(length(out1),1);
xrand=1:1:15;
for irat=1:length(yrandbyfun)
        N=100;m=out1(irat).parti_y{1}(1);
        yrandbyfun{irat}=randfun(xrand,m,N);
end
xrandmean=1:1:x(end);N=100;m=y(1);
yrandmean=randfun(xrandmean,m,N);
     
% -----by sampling
ysrandall=cellfun(@(x) mean(cell2mat(x),2), {out1.parti_y_rand},'UniformOutput',false);
% figure;imagesc(cell2mat(ysrand{1}))
ysrandmean=mean(cell2mat(ysrandall),2);
xsrand=out(1).trackx(:);
% 
figure('color','w');
plot([0;x],[0;y],'b.-');hold on;
plot([0;xsrand],[0;ysrandmean],'r.-');
plot([0;xrand(:)],[0;yrandmean(:)],'g.-');
legend({'Place cell recruit','Random Rate','Random'},'location','best','box','off')
xlabel('#Track');ylabel('%Cell');
% individual  animal
method={'hyperbolic'};
fitHyper=struct;
for irat=1:length(yall)
 tmp=eval_capacity_formula(yrandbyfun{irat}/100,xrand,method,0,{'#Track','%Cell'});
 tmp.datatype='Random'; if irat==1; fitHyper=tmp;else; fitHyper=cat(1,fitHyper,tmp);end
 tmp=eval_capacity_formula(ysrandall{irat}/100,xsrand,method,0,{'#Track','%Cell'});
 tmp.datatype='RandomRate'; fitHyper=cat(1,fitHyper,tmp);
 tmp=eval_capacity_formula(yall{irat}/100,x,method,0,{'#Track','%Cell'});
 tmp.datatype='Data'; fitHyper=cat(1,fitHyper,tmp);
end
% figure;plot(xsrand,ysrand,'k:');
% fitting
 method={'linear','log','log-add','power','exp','hyperbolic'};
 [fitres_rand0,test_residual_rand0,f0]=eval_capacity_formula(yrandmean/100,xrand,method,1,{'#Track','%Cell'});
 [fitres_rand,test_residual_rand,f1]=eval_capacity_formula(ysrandmean/100,xsrand,method,1,{'#Track','%Cell'});
 [fitres,test_residual,f2]=eval_capacity_formula(y/100,x,method,1,{'#Track','%Cell'});
 
%  [fitres_rand0,test_residual_rand0,f0]=eval_capacity_formula([0;yrandmean(:)/100],[0;xrand(:)],method,1,{'#Track','%Cell'});
%  [fitres_rand,test_residual_rand,f1]=eval_capacity_formula([0;ysrandmean/100],[0;xsrand],method,1,{'#Track','%Cell'});
%  [fitres,test_residual,f2]=eval_capacity_formula([0;y/100],[0;x],method,1,{'#Track','%Cell'});
 
 
%  savepdf(f0,fullfile(savepath,'recruitCurveRandomSampling_compareFiting.pdf'));
%  savepdf(f1,fullfile(savepath,'recruitCurveRandomSamplingByRate_compareFiting.pdf'));
%  savepdf(f2,fullfile(savepath,'recruitCurve_compareFiting.pdf'));
 
% savepdf(f3,fullfile(savepath,'recruitCurve_ncell_eachTrack.pdf'));
%% fig s5b ; random vs data
% methodshow={'linear','log','power','exponential','lomax','hyperbolic'};
methodshow={'linear','log','power','exponential','hyperbolic'};

f2=figure('color','w','position',[ 103    99   590   725]);
lcolor=jet(length(methodshow));
subplot(231)
plot(x,y,'ko');hold on;
for i=1:length(methodshow)
   fit1=fitres(strcmp({fitres.method},methodshow{i})); 
   plot(fit1.xi,fit1.yi*100,'color',lcolor(i,:));
end
set(gca,'ylim',[30,110],'xlim',[-1,16],'tickdir','out','box','off');
xlabel('#Track');ylabel('%Cell');title('Place cell recruitment');
subplot(232)
plot(xsrand,ysrandmean,'ko');hold on;h=[];
for i=1:length(methodshow)
   fit1=fitres_rand(strcmp({fitres.method},methodshow{i})); 
   h1=plot(fit1.xi,fit1.yi*100,'color',lcolor(i,:));
   h=[h;h1];
end
legend(h,methodshow,'location','best','box','off');
set(gca,'ylim',[30,110],'xlim',[-1,16],'tickdir','out','box','off');
xlabel('#Track');ylabel('%Cell');title('Random sampling');
subplot(233)
plot(xrand,yrandmean,'ko');hold on;h=[];
for i=1:length(methodshow)
   fit1=fitres_rand0(strcmp({fitres.method},methodshow{i})); 
   h1=plot(fit1.xi,fit1.yi*100,'color',lcolor(i,:));
   h=[h;h1];
end
set(gca,'ylim',[30,110],'xlim',[-1,16],'tickdir','out','box','off');
xlabel('#Track');ylabel('%Cell');title('Random');
subplot(212)
residualmse=zeros(length(methodshow),2);
residualmse_rand=zeros(length(methodshow),2);
residualmse_rand0=zeros(length(methodshow),2);
for i=1:length(methodshow)
    residualmse(i,1)=fitres(strcmp({fitres.method},methodshow{i})).residual_mean;
    residualmse(i,2)=fitres(strcmp({fitres.method},methodshow{i})).residual_se;
    residualmse_rand(i,1)=fitres_rand(strcmp({fitres_rand.method},methodshow{i})).residual_mean;
    residualmse_rand(i,2)=fitres_rand(strcmp({fitres_rand.method},methodshow{i})).residual_se;
    residualmse_rand0(i,1)=fitres_rand0(strcmp({fitres_rand0.method},methodshow{i})).residual_mean;
    residualmse_rand0(i,2)=fitres_rand0(strcmp({fitres_rand0.method},methodshow{i})).residual_se;
end
h1=errorbar(residualmse(:,1),residualmse(:,2));hold on;
h2=errorbar(residualmse_rand(:,1),residualmse_rand(:,2));
h3=errorbar(residualmse_rand0(:,1),residualmse_rand0(:,2));
set([h1,h2,h3],'linewidth',2)
xstr={'linear','log','power','poisson','hill'};
legend([h1,h2,h3],{'Place cell recruit','Random sampling','Random'},'location','best','box','off')
set(gca,'xlim',[0,length(methodshow)+1],'tickdir','out','box','off',...
    'xtick',1:length(methodshow),'xticklabel',xstr,'XTickLabelRotation',30);
ylabel('RSS');
setaxisformalAll

% hyperbolic parameter comparison
fit1=fitres(strcmp({fitres.method},'hyperbolic'));
fitr=fitres_rand(strcmp({fitres.method},'hyperbolic'));
fitr0=fitres_rand0(strcmp({fitres.method},'hyperbolic'));
c95=[fit1.fit95fun(fit1.fitparam,.95),fitr.fit95fun(fitr.fitparam,.95),fitr0.fit95fun(fitr0.fitparam,.95)];
f3=figure('color','w');
% subplot(211)
% plot(x,y,'bo');hold on;
% plot(xsrand,ysrandmean,'ro');
% h1=plot(fit1.xi,fit1.yi,'b');
% h2=plot(fitr.xi,fitr.yi,'r');
% legend([h1,h2],{'Place cell recruit','Random sampling'},'location','best','box','off')
% xlabel('#Track');ylabel('%Cell');
paramName={'n','C_{50}','R_{max}','C_{95}'};
for i=1:length(fit1.fitparam)
subplot(2,4,4+i);
plot([fit1.fitparam(i),fitr.fitparam(i),fitr0.fitparam(i)],'o-');
set(gca,'xtick',[1:3],'xlim',[0.5,3.5],'xticklabel',{'Data','RandomRate','Random'},...
    'box','off','xticklabelrotation',30);
ylabel(paramName{i});
end
i=4;
subplot(2,4,4+i);
plot(c95,'o-');
set(gca,'xtick',[1:3],'xlim',[0.5,3.5],'xticklabel',{'Data','RandomRate','Random'},...
    'box','off','xticklabelrotation',30);
ylabel(paramName{i});

setaxisformalAll

f4=figure('color','w');
subplot(221)
plot(x,y,'bo');hold on;
plot(xsrand,ysrandmean,'ro');
h1=plot(fit1.xi,fit1.yi,'b');
h2=plot(fitr.xi,fitr.yi,'r');
legend([h1,h2],{'Place cell recruit','Random sampling'},'location','best','box','off')
xlabel('#Track');ylabel('%Cell');
subplot(222)
plot(x,abs(y-ysrandmean),'ko');hold on;
xlabel('#Track');ylabel('|\Delta%Cell|');
subplot(223)
[tmp1,tmp2]=mean_se(abs(y-ysrandmean));
errorbarformal(1,tmp1,tmp2);
ylabel('|\Delta%Cell|');
setaxisformalAll
% savepdf(f3,fullfile(savepath,'recruitCurveDataVSRandom_compareFiting.pdf'));
% savepdf(f4,fullfile(savepath,'recruitCurveDataVSRandom_compareHRatioParameter.pdf'));
%   savepdf(f2,fullfile('/media/Big1/15tracks/figures/','recruitCurveDataVSRandom_compareModel.pdf'));
%   savepdf(f3,fullfile('/media/Big1/15tracks/figures/','recruitCurveDataVSRandom_compareHRatioParameter.pdf'));

% %  firingrate on fitting parameters
paramName={'n','C_{50}','R_{max}'};
rate_fitpa=[cat(1,out.ipath),cat(1,out.peakRateThr),cat(1, out.fitpa)];
peakRate=unique([out.peakRateThr]);
figure('color','w','position',[328    83   957   297]);
for i=1:3 
subplot(1,3,i)
y=[];
for j=1:5
ydata=rate_fitpa(rate_fitpa(:,1)==j,2+i);
y=[y;ydata(:)'];
end
[m,se]=mean_se(y);
shadedErrorBar(1:length(m),m,se,'b');
set(gca,'xtick',1:length(peakRate),'xticklabel',num2strcell(peakRate),'xlim',...
    [0,length(peakRate)+1],'tickdir','out','box','off')
xlabel('Peak Rate Threshold(Hz)');ylabel(paramName{i});
end
setaxisformalAll
% % % % % % savepdf(gcf,fullfile(savepath,'recruitCurve_PeakRate_HRatioParameter.pdf'));
%  savepdf(gcf,fullfile('/media/Big1/15tracks/figures/','recruitCurve_PeakRate_HRatioParameter.pdf'));
%%  recruitCurveDataVSRandom
% out1=out([out.peakRateThr]==2);
% yall=cellfun(@(x) x{1}(:),{out1.parti_y},'UniformOutput',false);
% y=mean(cell2mat(yall),2);
% x=out(1).trackx(:);
% 
% figure('color','w','position',[472    84   626   615]);
% subplot(121)
%  linec=jet(length(out1));
%     for i=1:length(out1)
% %        h=errorbar(out1(i).trackx,out1(i).parti_y{1},out1(i).parti_y{2},...
% %            'linewidth',1,'linestyle','-');hold on;
%        h=shadedErrorBar(out1(i).trackx,out1(i).parti_y{1},out1(i).parti_y{2},{'color',linec(i,:)},0.5);hold on;
% %        xrand=1:.1:15;N=100;m=out1(i).parti_y{1}(1);
% % %        yrand=(1-(1-m/N).^xrand)*100;
% %         yrand=randfun(xrand,m,N);
% %        plot(xrand,yrandbyfun{i},':','color',get(h,'color'))
% %        plot(xsrand,ysrandall{i},'-.','color',get(h,'color'))
%     end
%     xlabel('#Tracks');ylabel('%Cells');
%     subplot(122)
%     linec=hsv(3);
%      h=plot(x,y,'.','color',linec(1,:));hold on;
%      plot(fitres(strcmp({fitres.method},'hyperbolic')).xi,fitres(strcmp({fitres.method},'hyperbolic')).yi*100,...
%          'color',linec(1,:));
%      hr=plot(xrandmean,yrandmean,'.','color',linec(2,:));
%      plot(fitres_rand0(strcmp({fitres_rand0.method},'hyperbolic')).xi,...
%          fitres_rand0(strcmp({fitres_rand0.method},'hyperbolic')).yi*100,...
%          'color',linec(2,:));
%      hr2=plot(xsrand,ysrandmean,'.','color',linec(3,:));
%      plot(fitres_rand(strcmp({fitres_rand.method},'hyperbolic')).xi,...
%          fitres_rand(strcmp({fitres_rand.method},'hyperbolic')).yi*100,...
%          'color',linec(3,:));
%     xlabel('#Tracks');ylabel('%Cells');
%     legend([h;hr;hr2],{'place cell recruitment','random sampling','random sampling by rate'},'location','best');
%     setaxisformalAll
    
    
% savepdf(gcf,fullfile(savepath,'recruitCurveDataVSRandom.pdf'));
    
    
%     savepdf(gcf,fullfile('/media/Big1/15tracks/figures/','recruitCurveDataVSRandom.pdf'));
%     
%     close all
%     save(fullfile(savepath,'recruitCurveDataVSRandom.mat'));

% edit s_placecell_ntrack_participation_newrand

%% fig 5b 5c s5c place cell recruitment 
% edit s_placefield_continuous_run
%

load('/media/Big1/15tracks/recruitment/ContinuousRun.mat','out');



data=cat(1,out.data);
% fit poiss and gamma poiss to the distribution of  #plfield per cell
 [fitnegbino,cinegbino]=mle(data,'distribution','NegativeBinomial');
    [fitpoiss,cipoiss]=mle(data,'distribution','Poisson');
    
    a=fitnegbino(1);
    b=fitnegbino(2)/(1-fitnegbino(2));
    
    nplfx=[min(data):max(data)]';
    nplfy=histcounts(data,[nplfx(:)-.5;nplfx(end)+.5])';
    
    xi=nplfx(1):nplfx(end);
    yinb=pdf('NegativeBinomial',xi,fitnegbino(1),fitnegbino(2))*sum(nplfy);
    yips=pdf('Poisson',xi,fitpoiss(1))*sum(nplfy);
    
    % -------- recruitment 
    % data
     loc=[];
     for i=1:length(out);
         loc=[loc;
             out(i).ipath*ones(size(out(i).runseq,1),1),out(i).runseq(:,[1,2])]; 
         silentcell=out(i).silentcell;
         loc=[loc;
             out(i).ipath*ones(size(silentcell)),silentcell,Inf(size(silentcell))];
     end
     loc(:,end)=loc(:,end)-min(loc(:,end))+1;% start from 0
     locx=unique(loc(:,end));
     cumN=zeros(size(locx));
     for i=1:length(cumN)
        cumN(i)=length(unique(loc(loc(:,end)<=locx(i),[1,2]),'rows')); 
     end
    %  function        
    lomaxfun=@(x,a,b) (a.*(b.^a)./(x+b).^(a+1));
    lomaxcdf=@(x,a,b) 1-(b./(x+b)).^a;
    xicum=0:.001:1;
    y=lomaxfun(xicum,a,b);
%     yicum=cumsum(y)*mean(diff(xicum)); % literally compute cumulative df
    yicum=lomaxcdf(xicum,a,b);         % use formula to compute cdf
%     figure;plot(xicum,y,'.-');hold on;plot(xicum,yicum,'.-');
    ypoiss=fitpoiss*exp(-fitpoiss*xicum);
    yicumpoiss=cumsum(ypoiss)*mean(diff(xicum));
    
    figure('color','w','position',[205    69   671   390]);
    subplot(121)
    bar(nplfx,nplfy,'k');hold on;
    h1=plot(xi,yinb,'g');
    h2=plot(xi,yips,'r');
    set([h1,h2],'linewidth',2);
    set(gca,'xlim',[-2,max(nplfx)+2])
    legend([h1,h2],{'gamma poiss','poiss'},'location','best');
    xlabel('# place field per cell');ylabel('# cells');
    subplot(122)
    plot(locx,cumN/cumN(end),'k','linewidth',2);hold on;
     h1=plot(xicum*locx(end-1),yicum,'g');
     h2=plot(xicum*locx(end-1),yicumpoiss,'r');
%      legend([h1,h2],{'gamma poiss','poiss'},'location','best');
     set([h1,h2],'linewidth',1);
    set(gca,'xlim',[0,locx(end-1)])
    xlabel('Location(cm)');ylabel('Proportion of recruited cells');
    setaxisformalAll
%     savepdf(gcf,'/media/Big1/15tracks/recruitment/PlFieldsGammapoissDistribut.pdf');
%     savepdf(gcf,'/media/Big1/15tracks/figures//PlFieldsGammapoissDistribut.pdf');
    
    %----------------------------- connect all 15 tracks
    
    data=cat(1,out.data15);
    % fit poiss and gamma poiss to the distribution of  #plfield per cell
    [fitnegbino,cinegbino]=mle(data,'distribution','NegativeBinomial');
    [fitpoiss,cipoiss]=mle(data,'distribution','Poisson');
    
    negbinof=fit_sample_distribution(data,{'NegativeBinomial'});
   SSE=negbinof.residual; N=length(data);
sigmaSquare=SSE/N; % biased estimator of variance by ML
lnLg=-(N/2)*log(2*pi)-(N/2)*log(sigmaSquare)-1/(2*sigmaSquare)*SSE;
    poissf=fit_sample_distribution(data,{'Poisson'});
   SSE=poissf.residual; N=length(data);
sigmaSquare=SSE/N; % biased estimator of variance by ML
lnLp=-(N/2)*log(2*pi)-(N/2)*log(sigmaSquare)-1/(2*sigmaSquare)*SSE;
LR = 2*(lnLg - lnLp); % has a X2 distribution with a df equals to number of constrained parameters, here: 1
pval = 1 - chi2cdf(LR, 1);


    a=fitnegbino(1);
    b=fitnegbino(2)/(1-fitnegbino(2));
    
    nplfx=[min(data):max(data)]';
    nplfy=histcounts(data,[nplfx(:)-.5;nplfx(end)+.5])';
    
    xi=nplfx(1):nplfx(end);
    yinb=pdf('NegativeBinomial',xi,fitnegbino(1),fitnegbino(2))*sum(nplfy);
    yips=pdf('Poisson',xi,fitpoiss(1))*sum(nplfy);
    
    % -------- recruitment
    % data
    locx15=cell(length(out),1);
    locdata=cell(length(out),1);
    cumNp15=cell(length(out),1);
    for irat=1:length(out)
        silentcell=out(irat).silentcell15;
        loc=[out(irat).runseq15(:,[1,2]) ;
            silentcell,Inf(size(silentcell))];
        loc(:,end)=loc(:,end)-min(loc(:,end))+1;% start from 0
        locx=unique(loc(:,end));locx(~isfinite(locx))=[];
        cumN=zeros(size(locx));
        for i=1:length(cumN)
            cumN(i)=size(unique(loc(loc(:,end)<=locx(i),1),'rows'),1);
        end
        locdata{irat}=splitapply(@min,loc(:,end),findgroups( loc(:,1)));
        locx15{irat}=locx;
        cumNp15{irat}=cumN/out(irat).ncell_pyr;
    end
    locx=unique(cell2mat(locx15)')';
    locx(locx>min(cellfun(@max,locx15)))=[];
    cumNp=zeros(size(locx));
    for irat=1:length(cumNp15)
        tmp=interp1(locx15{irat},cumNp15{irat},locx);
        cumNp=cumNp+tmp;
    end
    cumNp=cumNp/length(cumNp15);
    
    %  function
    lomaxfun=@(x,a,b) (a.*(b.^a)./(x+b).^(a+1));
    lomaxcdf=@(x,a,b) 1-(b./(x+b)).^a;
    xicum=0:.001:1;
    y=lomaxfun(xicum,a,b);
    %     yicum=cumsum(y)*mean(diff(xicum)); % literally compute cumulative df
    yicum=lomaxcdf(xicum,a,b);         % use formula to compute cdf
    %     figure;plot(xicum,y,'.-');hold on;plot(xicum,yicum,'.-');
    ypoiss=fitpoiss*exp(-fitpoiss*xicum);
    yicumpoiss=cumsum(ypoiss)*mean(diff(xicum));
    
    yicum1=lomaxcdf(linspace(0,1,length(locx)),a,b)';
    locdata1=cell2mat(locdata);locdata1=locdata1(isfinite(locdata1));
    locdata1=locdata1(locdata1>=min(locx) & locdata1<=max(locx));
     [~,kspval,ksstat]=kstest(locdata1,'CDF',[locx,yicum1/yicum1(end)])
%    
%     figure;plot(locx,yicum1);hold on;plot(locx,cumNp);cdfplot(locdata1);
%     locdata1=cell2mat(locdata);locdata1=locdata1(isfinite(locdata1));
%     [fitnegbino1,cinegbino1]=mle(locdata1,'distribution','NegativeBinomial');
%     test_cdf=makedist('NegativeBinomial','R',fitnegbino1(1),'p',fitnegbino1(2));
%     [~,kspval,ksstat]=kstest(locdata1,'CDF',test_cdf)
    
    figure('color','w','position',[205    69   671   390]);
    subplot(121)
    nplfy_norm=nplfy/sum(nplfy);
    bar(nplfx,nplfy_norm,'k');hold on;
    h1=plot(xi,yinb/sum(yinb),'g');
    h2=plot(xi,yips/sum(yips),'r');
    set([h1,h2],'linewidth',2);
    set(gca,'xlim',[-2,max(nplfx)+2])
    legend([h1,h2],{'gamma poiss','poiss'},'location','best');
    xlabel('# place field per cell');ylabel('# cells');
    subplot(122)
    plot(locx,cumNp,'k','linewidth',2);hold on;
    h1=plot(xicum*locx(end),yicum,'g');
    h2=plot(xicum*locx(end),yicumpoiss,'r');
    %      legend([h1,h2],{'gamma poiss','poiss'},'location','best');
    set([h1,h2],'linewidth',1);
    set(gca,'xlim',[0,locx(end-1)],'ylim',[0,1])
    xlabel('Location(cm)');ylabel('Proportion of recruited cells');
    setaxisformalAll
%     savepdf(gcf,'/media/Big1/15tracks/recruitment/PlFieldsGammapoissDistribution_All15.pdf');
    
%      savepdf(gcf,'/media/Big1/15tracks/figures/PlFieldsGammapoissDistribution_All15.pdf');
%% fig 5a; place cell recruitment 
clear
 load('/media/Big1/15tracks/recruitment/recruitCurveDataVSRandom.mat');
 
 
out1=out([out.peakRateThr]==2);
yall=cellfun(@(x) x{1}(:),{out1.parti_y},'UniformOutput',false);
y=mean(cell2mat(yall),2);
x=out(1).trackx(:);

f1=figure('color','w','position',[472    84   626   615]);
subplot(121)
 linec=jet(length(out1));
    for i=1:length(out1)
%        h=errorbar(out1(i).trackx,out1(i).parti_y{1},out1(i).parti_y{2},...
%            'linewidth',1,'linestyle','-');hold on;
       h=shadedErrorBar(out1(i).trackx,out1(i).parti_y{1},out1(i).parti_y{2},{'color',linec(i,:)},0.5);hold on;
%        xrand=1:.1:15;N=100;m=out1(i).parti_y{1}(1);
% %        yrand=(1-(1-m/N).^xrand)*100;
%         yrand=randfun(xrand,m,N);
%        plot(xrand,yrandbyfun{i},':','color',get(h,'color'))
%        plot(xsrand,ysrandall{i},'-.','color',get(h,'color'))
    end
    xlabel('#Tracks');ylabel('%Cells');
    subplot(122)
    linec=hsv(3);
     h=plot(x,y,'.','color',linec(1,:));hold on;
     plot(fitres(strcmp({fitres.method},'hyperbolic')).xi,fitres(strcmp({fitres.method},'hyperbolic')).yi*100,...
         'color',linec(1,:));
     hr=plot(xrandmean,yrandmean,'.','color',linec(2,:));
     plot(fitres_rand0(strcmp({fitres_rand0.method},'hyperbolic')).xi,...
         fitres_rand0(strcmp({fitres_rand0.method},'hyperbolic')).yi*100,...
         'color',linec(2,:));
     hr2=plot(xsrand,ysrandmean,'.','color',linec(3,:));
     plot(fitres_rand(strcmp({fitres_rand.method},'hyperbolic')).xi,...
         fitres_rand(strcmp({fitres_rand.method},'hyperbolic')).yi*100,...
         'color',linec(3,:));
    xlabel('#Tracks');ylabel('%Cells');
    legend([h;hr;hr2],{'place cell recruitment','random sampling','random sampling by rate'},'location','best');
    setaxisformalAll
 
 
 
% hyperbolic parameter comparison
fit1=fitres(strcmp({fitres.method},'hyperbolic'));
fitr=fitres_rand(strcmp({fitres_rand.method},'hyperbolic'));
fitr0=fitres_rand0(strcmp({fitres_rand0.method},'hyperbolic'));
c95=[fit1.fit95fun(fit1.fitparam,.95),fitr.fit95fun(fitr.fitparam,.95),fitr0.fit95fun(fitr0.fitparam,.95)];
f3=figure('color','w');
paramName={'n','C_{50}','R_{max}','C_{95}'};
for i=1:length(fit1.fitparam)
subplot(2,2,i);
plot([fit1.fitparam(i),fitr.fitparam(i),fitr0.fitparam(i)],'o-');
set(gca,'xtick',[1,2,3],'xlim',[0.5,3.5],'xticklabel',{'pc recruit','Random by Rate','Random'},...
    'box','off','xticklabelrotation',30);
ylabel(paramName{i});
end
i=4;
subplot(2,2,i);
plot(c95,'o-');
set(gca,'xtick',[1,2,3],'xlim',[0.5,3.5],'xticklabel',{'pc recruit','Random by Rate','Random'},...
    'box','off','xticklabelrotation',30);
ylabel(paramName{i});

setaxisformalAll


n951=floor([fitHyper(strcmp({fitHyper.datatype},'Data')).ntrackfit95]);
n95r=floor([fitHyper(strcmp({fitHyper.datatype},'Random')).ntrackfit95]);
signrank(n951,n95r)

%   savepdf(f1,fullfile('/media/Big1/15tracks/figures/','recruitCurveDataVSRandom.pdf'));
%     
%   savepdf(f3,fullfile('/media/Big1/15tracks/figures/','recruitCurveDataVSRandom_compareHRatioParameter.pdf'));
  
  %% fig 5d 5e  s5h s5i  ; placeCellParticipation sleepFiringRate
  
clear
load('/media/Big1/15tracks/SleepFiringRate/firingrateWithBurst_population170220131512.mat')


mse=cell(4,3);
for i=1:3
    for isleep=1:4
        xy=data4sleep{isleep}(:,[i,4]);
        [ind,ntrx]=findgroups(xy(:,2));
        [m,se]=splitapply(@mean_se,xy(:,1),ind);
        mse{isleep,i}=[ntrx(:),m(:),se(:)];
    end
end

[parcc,parp]=partialcorr(data4sleep{1}(:,4),data4sleep{4}(:,1),data4sleep{1}(:,1));
diffsleep14=data4sleep{4}(:,1)-data4sleep{1}(:,1);
[ind,ntrx]=findgroups(data4sleep{1}(:,4));
[diffsleep14m,diffsleep14se]=splitapply(@mean_se,diffsleep14,ind);


f1=figure('color','w','position',[205   228   829   724]);
for i=1:3
for isleep=1:4
    subplot(3,4,isleep+(i-1)*4);
    plot(data4sleep{isleep}(:,4),data4sleep{isleep}(:,i),'r.');hold on
    shadedErrorBar(mse{isleep,i}(:,1),mse{isleep,i}(:,2),mse{isleep,i}(:,3),{},0.5);
    xlabel(dataStr{4});ylabel(dataStr{i});
    title(sprintf('Sleep%d\nr=%.2g,p=%.2g',isleep,ccrp_ntrack(isleep,1,i),ccrp_ntrack(isleep,2,i)));
    lsline;
end
end
setaxisformalAll(10)



f2=figure('color','w','position',[  117   257   844   710]);
ecolor=hsv(4);
for i=1:3 
subplot(3,3,i)
   shadedErrorBar(1:4, data4sleepmse(:,i,1),data4sleepmse(:,i,2),{'color','r'});
xlabel('Sleep');ylabel(dataStr{i});set(gca,'xlim',[0,5],'xtick',1:4);

subplot(3,3,3+i)
h=[];
for isleep=1:4
xy=data4sleep{isleep}(:,[i,4]);
[ind,ntrx]=findgroups(xy(:,2));
[m,se]=splitapply(@mean_se,xy(:,1),ind);
h1=shadedErrorBar(ntrx,m,se,{'color',ecolor(isleep,:)},0.5);hold on;h=[h;h1.mainLine];
end
xlabel('# of tracks');ylabel(dataStr{i});
set(gca,'xlim',[0,16],'xtick',[0,5:5:20]);
if i==1
legend(h,{'Sleep1';'Sleep2';'Sleep3';'Sleep4'},'location','best','box','off');
end
subplot(3,3,6+i)
plot(1:4, ccrp_ntrack(:,1,i),'o-');
xlabel('Sleep');ylabel(sprintf(['Correlation\n(',dataStr{i},'~Ntrack)']));
set(gca,'xlim',[0,5],'xtick',1:4);
end
setaxisformalAll(10)

f3=figure('color','w');
subplot(121); i=1;
h=[]; m0=zeros(16,4);se0=zeros(16,4);
for isleep=[1:4]
xy=data4sleep{isleep}(:,[i,4]);
[ind,ntrx]=findgroups(xy(:,2));
[m,se]=splitapply(@mean_se,xy(:,1),ind);
m0(:,isleep)=m; 
se0(:,isleep)=se;
h1=shadedErrorBar(ntrx,m,se,{'color',ecolor(isleep,:)},0.5);hold on;h=[h;h1.mainLine];
end
xlabel('# of tracks');ylabel(dataStr{i});
set(gca,'xlim',[0,16],'xtick',[0,5:5:20]);
legend(h,{'Sleep1';'Sleep2';'Sleep3';'Sleep4'},'location','best','box','off');
subplot(122);
% [parcc,parp]=partialcorr([0:15]',m0(:,4),m0(:,1));
% diffsleep14=m0(:,4)-m0(:,1);
shadedErrorBar(0:15,diffsleep14m,diffsleep14se);
title(sprintf('Partial corr R=%.2g,P=%.2g',parcc,parp));
xlabel('# of tracks');ylabel(['\Delta ',dataStr{i},' S4-S1']);
setaxisformalAll(10)
% saveallfig2pdf([savepath,'/FiringBurstSleep_vsNumberOfTrackParticipation.pdf']);

%   savepdf(f1,fullfile('/media/Big1/15tracks/figures/','placeCellParticipation_sleepFiringRate_scatter.pdf'));
%   savepdf(f2,fullfile('/media/Big1/15tracks/figures/','placeCellParticipation_sleepFiringRate_summary.pdf'));
%     
%   savepdf(f3,fullfile('/media/Big1/15tracks/figures/','placeCellParticipation_sleepFiringRate_sleep14Diff.pdf'));
%% fig s5f; frequency of each longest common sequence

load('/media/Big1/15tracks/LongestCommonSeq/lcsoutput.mat');
good=struct;
for irat=1:length(out)
    
    occur=cell(size(out(irat).lcs1tr,1),4); % [1 cellid, 2run pair id; 3 number of pair; 4 run pairs; 5 tracks]
    for i1=1:size(out(irat).lcs1tr,1)
        occur{i1,1}=out(irat).lcs1tr{i1,3};
        tmp=[];
        for i2=1:size(out(irat).lcs1tr,1)
            if contains(char(out(irat).lcs1tr{i2,3}), char(out(irat).lcs1tr{i1,3}))
                
                tmp=[tmp;i2];
            end
        end
        occur{i1,2}=tmp;
        occur{i1,3}=length(tmp);
        occur{i1,4}=out(irat).pair1tr_tr(tmp,:);
        
        tmp1=out(irat).pair1tr_tr(tmp,:);
        
        tmp1=[[tmp1(:,1),tmp1(:,3)];tmp1(:,[2,3])];
        occur{i1,5}=unique(tmp1,'rows');
        
    end
    
    goodlcsind=find(cell2mat(occur(:,3))>0 & out(irat).lcs_len >=2  & out(irat).pval_lcs1tr<0.05);
    
    goodlcs=occur(goodlcsind,:);
    
    goodlcsnote=['1 cellid, 2run pair id; 3 number of pair; 4 run pairs; 5 tracks'];
    
    good(irat).goodlcs=goodlcs;
    good(irat).goodlcsnote=goodlcsnote;
    good(irat).ratname=out(irat).ratname;
end

goodlcs=cat(1,good.goodlcs);

goodlen=cellfun(@length,goodlcs(:,1));
goodntrack=cellfun(@(x) length(unique(x(:,1))),goodlcs(:,5));

[n2,edge1,edge2]=histcounts2(goodlen,goodntrack);

x1=edge2x(edge1);
x2=edge2x(edge2);
f1=figure('position',[   366   557   874   421]);
subplot(121)
imagesc(x2,x1,n2);
xlabel('# of tracks');ylabel('length of sequence');
title('# of longest common sequence');
colormap copper
subplot(122)
h=histogram(goodlen,0:max(goodlen)*1.5);
set(h,'facecolor','k','facealpha',1);
xlabel('Length of LCS');
ylabel('Count of LCS');
setaxisformalAll

[m,se]=mean_se(goodlen)

savepdf(f1,'/media/Big1/15tracks/figures/LCS_length.pdf')
%% longest common seq vs orientation

load('/media/Big1/15tracks/LongestCommonSeq/lcsoutput_orientation.mat');

f1=figure('position',[ 258    85   647   630]);
subplot(2,2,1)
errorbarformal(1:4,squeeze(lcsdata_meanse(1,1,:)),squeeze(lcsdata_meanse(1,2,:)));
sigstar({[1,2],[3,4]},[pval_ratio_headdir,pval_ratio_headori]);
set(gca,'xtick',1:4,'xticklabel',ind_note,'XTickLabelRotation',30);
ylabel('length ratio');
subplot(2,2,2)
errorbarformal(1:4,squeeze(lcsdata_meanse(2,1,:)),squeeze(lcsdata_meanse(2,2,:)));
sigstar({[1,2],[3,4]},[pval_len_headdir,pval_len_headori]);
set(gca,'xtick',1:4,'xticklabel',ind_note,'XTickLabelRotation',30);
ylabel('length ratio');

subplot(2,2,3);
plot(lcs_count([1,2],:),'o-');
set(gca,'xlim',[0,3],'xtick',1:2,'xticklabel',ind_note(1:2),'XTickLabelRotation',30);
ylabel('proportion of sig');
subplot(2,2,4);
plot(lcs_count([3,4],:),'o-');
set(gca,'xlim',[0,3],'xtick',1:2,'xticklabel',ind_note(3:4),'XTickLabelRotation',30);
legend('location','best');
ylabel('proportion of sig');
setaxisformalAll(12)
suptitle(['All rats']);

savepdf(f1,'/media/Big1/15tracks/figures/LCS_ori_dir.pdf')

%% fig s5g; longest common seq vs tuplet
load('/media/Big1/15tracks/LongestCommonSeq/tupletInLCs.mat','tuprand','good',...
    'tuplcs','tuplcsrand');

typestr={'Cells','Same order','Rev order'};
data={[tuplcs.ratio_tup_in],[tuprand.ratio_tup_in],[tuplcsrand.ratio_tup_in]};
[m,se,p]=mean_se_cell(data);
[~,gpair,p]=test_multigroup_pairwise(cell2mat(data')');
data={[tuplcs.ratio_tup_order1],[tuprand.ratio_tup_order1],[tuplcsrand.ratio_tup_order1]};
[m1,se1,p1]=mean_se_cell(data);
[~,gpair1,p1]=test_multigroup_pairwise(cell2mat(data')');
data={[tuplcs.ratio_tup_order2],[tuprand.ratio_tup_order2],[tuplcsrand.ratio_tup_order2]};
[m2,se2,p2]=mean_se_cell(data);
[~,gpair2,p2]=test_multigroup_pairwise(cell2mat(data')');

% f1=figure('position',[175         100        1333         367]);
% subplot(131);
% errorbarformal(1:length(m),m*100,se*100,{'data','randomSleep','non-com seq'});
% sigstar(gpair,p(:,end));
% title(typestr{1}); ylabel('Proportion of tuplets(%)');
% set(gca,'XTickLabelRotation',30);
% subplot(132);
% errorbarformal(1:length(m),m1*100,se1*100,{'data','randomSleep','non-com seq'});
% sigstar(gpair1,p1(:,end));
% title(typestr{2}); ylabel('Proportion of tuplets(%)');
% set(gca,'XTickLabelRotation',30);
% subplot(133);
% errorbarformal(1:length(m),m2*100,se2*100,{'data','randomSleep','non-com seq'});
% sigstar(gpair2,p2(:,end));
% title(typestr{3}); ylabel('Proportion of tuplets(%)');
% set(gca,'XTickLabelRotation',30);
% setaxisformalAll
m=m(1:2);m1=m1(1:2);m2=m2(1:2);
se=se(1:2);se1=se1(1:2);se2=se2(1:2);
f1=figure('position',[175         100        1333         367]);
subplot(131);
errorbarformal(1:length(m),m*100,se*100,{'data','randomSleep'});
sigstar(gpair{1},p(1,end));
title(typestr{1}); ylabel('Proportion of tuplets(%)');
set(gca,'XTickLabelRotation',30);
subplot(132);
data={[tuplcs.ratio_tup_order1],[tuprand.ratio_tup_order1]};
[ydata,gr]=cell2group(data);
plot_errorbar_with_data(ydata*100,gr)
% errorbarformal(1:length(m),m1*100,se1*100,{'data','randomSleep'});
sigstar(gpair1{1},p1(1,end));
title(typestr{2}); ylabel('Proportion of tuplets(%)');
set(gca,'XTickLabelRotation',30);
subplot(133);
errorbarformal(1:length(m),m2*100,se2*100,{'data','randomSleep'});
sigstar(gpair2{1},p2(1,end));
title(typestr{3}); ylabel('Proportion of tuplets(%)');
set(gca,'XTickLabelRotation',30);
setaxisformalAll

% savepdf(f1,'/media/Big1/15tracks/figures/LCS_tuplet.pdf')
savepdf(f1,'/media/Big1/15tracks/figures/LCS_tuplet_dot.pdf')